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' f-i ' Abstract 

Oh 

Q A robust second order, shock-capturing numerical scheme for multi-dimensional 

special relativistic magnetohydrodynamics on computational domains with adap- 
tive mesh refinement is presented. The base solver is a total variation diminishing 
Lax-Friedrichs scheme in a finite volume setting and is combined with a diffusive 
approach for controlling magnetic monopole errors. The consistency between the 
^ primitive and conservative variables is ensured at all limited reconstructions and 

the spatial part of the four velocity is used as a primitive variable. Demonstra- 
f- — five relativistic examples are shown to validate the implementation. We recover 

<^ known exact solutions to relativistic MHD Riemann problems, and simulate the 

t^^ shock-dominated long term evolution of Lorentz factor 7 vortical flows distorting 

^^ magnetic island chains. 
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1 Introduction 



Relativistic flows of magnetized plasmas are observed in a wide variety of 
astrophysical objects. At galactic scales, accretion of matter around black holes 
in Active Galactic Nuclei (AGN) in addition results in well collimated jets 
emitted along the rotation axis. To explain their observed superluminal plasma 
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motion at the parsec length scales, flows with Lorentz factors of more than 10 
are needed. The observed synchrotron emission indicates that these jets are 
pervaded by magnetic flelds. Even more powerful processes are at play in the 
highly relativistic blast waves associated with Gamma Ray Bursts (GRBs). 
Here, the plasma flows easily reach Lorentz factors of 100 or even higher. The 
morphology and time evolution of these and other astrophysical objects often 
involve strong shocks and complex magnetic fleld topologies. To compute this 
kind of challenging relativistic plasma dynamics, the combination of Adaptive 
Mesh Reflnement (AMR) with a robust, shock-capturing numerical method is 
therefore indispensible. 

With growing interest in relativistic astrophysical phenomena, various efforts 
are ongoing to develop numerical special relativistic hydrodynamic and mag- 
netohydrodynamic codes. Signiflcant progress was achieved in the last decade 
with the development of conservative shock-capturing schemes, which use ei- 
ther exact Riemann solvers or approximate Riemann solvers, or more robust 
central type schemes, in relativistic hydrodynamics |14|16|13|IT1lll|29|31j (for 
a contemporary review see [29j) and in relativistic magnetohydrodynamics 
PHH26|12|[7|27p32|33|20j . The study of relativistic hydrodynamic fluids cur- 
rently starts to beneflt also from using spatial and temporal adaptive tech- 
niques, or Adaptive Mesh Reflnement (AMR) [iill30lli3] . To date, AMR was 
incorporated only in some works, both in special relativistic magnetohydrody- 
namics [7| and in general relativistic magnetohydrodynamics pi^llTH] . At the 
same time, various authors, including [30j, have shown that AMR is impera- 
tive to simulate extreme astrophysical phenomena, such as those encountered 
in GRBs. 

In order to handle GRB and other extreme relativistic flow regimes, the em- 
ployed solver needs to be very robust under a wide variety of plasma condi- 
tions. Therefore, we decided to use a conservative discretization which does 
not fully exploit the detailed solution knowledge of the Riemann problem 
at each cell interface, deflned by two constant states in contact. This con- 
trasts with a true Godunov scheme, which would compute the flux across 
cell interfaces from the exact solution to the local Riemann problem. While 
fairly recently [I7j, an exact Riemann problem solver for the RMHD equa- 
tions became available, the nonlinear iteration involved would make it a fairly 
computationally costly ingredient for a multidimensional code. As pointed 
out further, we already face a similar Newton-Raphson procedure to deduce 
primitive from conservative variables in every grid point. In fact, we follow 
the trend away from using exact or approximate Riemann problem solvers 
in most codes in use to date. In the works cited above, Komissarov [21] al- 
ready used a linearized (approximate) Riemann solver, and even fell back on 
an HLLE variant [T9p5] which only uses the fastest wave speeds. Koldoba 
et al. [2S] and Balsara |1] also employ an approximate, linearized Riemann 
solver and presented numerical solutions for stringent ID test problems. All 



other works mentioned [T2|7|27II32||20] use central-type schemes, of either the 
local Lax-Friedrichs, HLL, or HLLC variant, suitably adopted to the relativis- 
tic MHD regime. We will use the simplest of these in our implementation, 
since in combination with automated grid adaptivity, we rely on convergence 
and robustness of the underlying discretization, while accuracy is most easily 
gained by raising the effective resolution employed. 

The macroscopic dynamics of astrophysical objects with relativistic plasma 
flows are governed by the conservation of particle and momentum-energy, 
along with the Maxwell equations to account for the involvement of the mag- 
netic field. In Sect. [2| we recast the RMHD equations in conservation form. 
Details on the implementation of the high-resolution shock-capturing scheme 
and the adaptive mesh refinement in our code can be found in Sects. [3]j4j The 
code is validated against several known ID analytical solutions in Sect. [5] and 
two multi-dimensional relativistic flow simulations are performed. In Sect. |6} 
a summary of the paper containing our main conclusions and an outlook for 
future applications of the code is given. 



2 Governing Equations 



We give a brief outline of the derivation of the relativistic magnetohydrody- 
namic (RMHD) equations and restrict the analysis to special relativity only. 
For a more extensive account we refer to [21]. Along the way, we point out 
minor differences in the choice of conservative variables, as well as a possibility 
to exploit an entropy related variable. The dynamics follow from the Maxwell 
equations and the conservation of particle number density and energy-stress. 
We will fix a Lorentzian reference frame to perform the computations. Due 
to Lorentz contraction the fluid volume elements are contracted so that the 
particle number density is multiplied by the Lorentz factor F = (1 — v^)~^/^. 
Here, v = {vcc,Vy,Vz) is the spatial velocity three-vector. We measure in the 
often used scaling for which the speed of light is unity c = 1. The particle 
number conservation equation can generally be written as 

da^ipU^ = 0, (1) 

where p is the proper mass density, i.e., the density in that Lorentz frame in 
which the fluid is locally at rest, u°' = F(l,v) is the four velocity that has 
to satisfy the constraint UaU"" = — 1, and Greek indices denote four vectors. 
We will use Latin indices to denote three-vectors. The Einstein summation 
convention for repeated indices is assumed. 

The antisymmetric dual field tensor J^'^^ can be defined in terms of the mag- 
netic field B and electric field E as measured in the Lorentzian reference frame 



by :r°/5 = -:r/3«, jF^^ = B\ P^ = -e'^^Ek, where e'^'' is the Levi-Cevita sym- 
bol. The homogeneous Maxwell equations can be expressed in terms of this 
dual tensor as 

a^jr°/3 = 0. (2) 



By subsequently assuming the fluid to be perfectly conducting, so that the 
CO moving electric field vanishes orE + vxB = 0, and by introducing the 
four-vector 

r = j^^Pup = (rv ■ B, B/r + rv ■ Bv), (3) 



the induction equation and V ■ B = transform to 

d^iu^'b'^ - b^u^) = 0. (4) 

Finally, conservation of energy and momentum follow from 

dp{Tf + T:i) = 0. (5) 

Here 



Tf = phu''u^+pr]''^, (6) 

T^^ = |6|2(m"m/3 + ^T]''^) - b'^bP, (7) 

are the fluid and electromagnetic energy-stress tensor, rf^ = diag(— 1, 1, 1, 1) 
is the metric tensor in fiat Minkowski space-time, and h is the relativistic 
specific enthalpy of the fluid in the local comoving frame. In the derivation of 
the electromagnetic tensor, use has been made of the Maxwell equations. 

To make the equations more useful for the numerical approach based on tem- 
porally advancing conservation laws, we split the equations for the mass den- 
sity ([I|, momentum and energy ([s]), and the magnetic field Q in the space 
and time coordinates: 



dt{Tp) + d,{Tpv')=0, (8) 

dt{T^ph,,,v^ - b'^V) + di{V^pK,,v'v^ + p,o,5'^ - UV) = 0, (9) 

dtiV^pK,, - j9tot - [b^f) + d,{T^pK,,v' - b%^) = 0, (10) 

dtB^ + di{v'B^ - B'v^) = 0, diB' = 0. (11) 

In these equations, we have introduced the relativistic magnetic and total 
pressure 



Pmag=^(^ + (V-B)2), (12) 

Ptot=P + Pma.g, (13) 

where p is the thermal pressure in the local comoving frame. The total rela- 
tivistic specific enthalpy is given by 

h,, = h+^-P^. (14) 



The equations (^8|-(11) are explicitly in conservation form. We actually com- 



bine the first and third equation as written above to arrive at 

dtV + V- F(U) = 0, (15) 

where 



V={D,S,E,Bf 



(rp, (e + i?2)v - (v ■ B)B, 



^+^ + l{v^B^-{^r.By)~p-D,B\ , (16) 

F = f Dv, (e + fi')vv - ?J - (v ■ B)(Bv + vB) + Jptot, 

^v + ptotV-(v-B)B,vB-Bv)^, (17) 

Here, TD = T'^p is the mass density in the inertial reference frame and 
^ = r'^ph is a measure for the enthalpy. The variable E has a contribution 
from the rest mass density subtracted from the total energy density as written 
for the laboratory frame, and represents a small difference between our vari- 
ables set and the one used in most implementations [T2I211I1]- This makes the 
current set of variables reduce to the usual set employed in non-relativistic 
MHD computations in the limit of F -^ 1. This also helps to avoid poten- 
tial numerical problems with negative pressures, in cases when the rest mass 
contribution is a dominant term in the energy density. In the case of an ideal 
equation of state with a constant polytropic index 7, the variable ^ becomes 



The quantity ^ (and ph) reduces to the non-relativistic enthalpy, except for 
the inclusion of the rest mass contribution. 



It is noted that instead of choosing the energy density £" as a conserved vari- 
able, we can also switch to DS, where S = p/p^ is the specific entropy. The 
energy equation is then replaced by 

dt{DS) + V ■ (DSy) = 0, (19) 



expressing the conservation of entropy for ideal RMHD. This conservation law 
also follows from UadpT^^ = hadpT""^ = 0, so that UadpT^ = Uadf^T"^ = 
leading after some algebra to the conservation of entropy. The energy equation 
has the distinct disadvantage that numerical calculations with small kinetic 
pressure compared to magnetic pressure can, depending on the V ■ B strategy 
in use, easily result in negative pressure. This problem could be circumvented 
by switching to the entropy equation. The use of this entropy conservation law 
was also mentioned in [20] , where it is advised for use in applications without 
shocks. Even on a fixed grid, a strategy for using this equation locally instead of 
the energy density law would need to be provided, explaining when positivity 
is preferred above strict conservation. Examples of such switching strategies 
can be found in challenging cosmological (non-relativistic) hydro codes, such 
as in [38|9] . In the shock-dominated simulations described below, we as yet 
did not need to use such a switching strategy. However, the grid adaptive code 
employed here allows to use a related idea, but only at the time of restriction 
and prolongation actions. Before every restriction/prolongation, we have the 
option to locally convert to the set of conserved variables [D, S, DS, B), keep- 
ing both conservation and positivity garantueed in these operations. After the 
coarsening or refining action, we then revert to the original conservative set 
employed in the time integration routines. This makes this strategy particular 
to AMR computations which can involve strong shocks, such as some of the 
tests presented below. While we did use this in some of our non-relativistic 
MHD computations so far, the tests mentioned below were performed without 
the need to invoke this entropy strategy. 



3 Relativistic Shock Capturing Scheme 

3.1 Switching from Conservative to Primitive Variables 



To advance the set of conservation equations (15-17) in time, we need to calcu- 
late the fiuxes. The latter are obtained from the primitive variables (p, v, p, B). 
While the determination of the primitive variables for non-relativistic MHD is 
a straightforward algebraic manipulation, the transformation of the variables 
for the relativistic MHD equations needs a root finding algorithm. The anal- 
ysis is facilitated by the introduction of the auxiliary variables (r, ^), which 



are to be determined from the conservative variables {D,S,E,B) only. Once 
the auxiliary variables are known, the construction of the primitive variables 
will be straightforward. 



We follow the method given in [7] to determine ^. From the definition of the 



conservative variables (16) and (F,^), we obtain the primitive variables 






(20) 



and 



P = D/T, 



P 



7 - 1 (^ - TD) 



7 



r2 



(21) 



Using the definition of E in Eq. (16) and the expression of ^ (18), it follows 
that ^ is the root of 



m)-^-'-^^^^-E-D,B^^l 



B^ (S-B)2 



0, (22) 



where the second term is the kinetic pressure p. The still unknown Lorentz 
factor r = r(S,B;,^) contains once more the variable S, as follows from Eq. 



(20) 



1 



1 - w^ = 1 - 



|S + rHS-B)B| 

(e + B^f 



(23) 



In the root finding algorithm, ^ has to be found as a zero of Eq. (22) with the 



help of the augmented equation ( 23 ) . 



The root ^ is found by means of a Newton-Raphson method. In our imple- 
mentation, the brackets for the Newton-Raphson method are found as follows: 
We constrain the pressure by a given lower bound p^, so that the values of ^ 
are restricted from below by 



e = r2(p + -^p) >D + -^p, = ^1, 



7-r 



7-r 



(24) 



since F > 1. The velocity that follows from this bound on ^i is constrained to 



a value below the speed of light, that means, by using Eq. (23) 



^^'(ei) < (1 - dv^y 



(25) 



where dv^ represents a given threshold. If ^i does not satisfy this condition, 
then find a new ^i via a Newton-Raphson procedure from t'^(^) = (1 — dv^Y- 
The obtained ^i is our lower bound. As a maximal bound for the bracketing 
we choose ^2 = E + D +p^. Next, we check the sign of /(G) and f{^2)- If they 
are the same, then the brackets are wrong, and a new guess for the brackets is 
found by successively replacing ^ by ^ and ^ by 2^. We follow this procedure 
till we have correct brackets. Finally, we apply the Newton-Raphson method 
to find C, and F, followed by a consistency check to ensure that v < 1 and 

P >Pe- 



3.2 Determining the Characteristic Speeds 



The shock-capturing scheme used in this work needs the (maximal) charac- 
teristic speeds of a given state of a fluid element. Since the equations are 
formulated in the inertial reference frame, we need to determine the charac- 
teristic wave speeds in this frame. 



For each spatial dimension i the RMHD equations (pj)-(ll) yields seven char- 
acteristic speeds. There is one entropy wave speed Ae corresponding to the 
passive advection of entropy disturbances in ideal RMHD. The other char- 
acteristic wave speeds, however, come in pairs, namely, two slow magneto- 
acoustic waves speed Af , two Alfven wave speeds A^, and two fast magneto- 
acoustic wave speeds Ap. The characteristic speeds are ordered according to 
the sequence of inequalities 

Ap < Aa < Ag < Ae < A^ < A^ < Aj^, (26) 



similar to the non-relativistic MHD equations. The entropy wave is just prop- 
agating with the fluid velocity Ae = f «. The Alfven disturbances propagate at 
speeds 

^A = ^.±4^=7=§^— ^, (27) 

r^ vp^tot ± (v ■ B) 



which includes relativistic corrections. The magneto-acoustic characteristic 
speeds follow from the quartic polynomial 



ph{\-l)T\X-v,Y-{l-\ 
ci 



2p 
C 



T'{ph+^-^){x-v,r- 



s 



F(vB)(A-^;,)-f^]%0, (2J 



and involves the relativistic hydro dynamic speed of sound Cs = J'yp/{ph). 
These magneto-acoustic wave speeds can be found algebraically, but the for- 
mulae are not numerically efficient and often susceptible to round-off errors. 
Instead we use Laguerre's method to find the four magneto-acoustic roots in 
the actual implementation. All roots are bound by the speed of light, so that 
they must lie in the interval ] — 1,1[. The roots can be located close to each 
other, especially when they are of the order unity. Making the transforma- 
tion /i = 1/(1 — A), we obtain a quartic polynomial for /i. The roots are now 
better separated and are in the interval ]0.5, cxo[. If the requested accuracy 
is not obtained via Laguerre's method, a root polishing is performed by the 
Newton-Raphson method. 

The above mentioned scheme for finding characteristic speeds can possibly be 
improved if the magneto-acoustic roots are almost degenerate. In that case 
root polishing by the Newton-Raphson procedure can unintentionally make 
the roots degenerate. It is then preferable to switch to Maehly's procedure 
(see [37]). Another improvement can be made by using the transformation 
/i = 1/(1 — \ + Vi) to reduce the possibility of excessively large roots. Another 
option which we implemented in our code is to compute the zeros in terms 
H = r(A — fj), which is somewhat better behaved for very large Lorentz factor 
flows. 



3.3 Total Variation Diminishing Lax-Friedrichs Scheme 



In the present paper, we employ the Total Variation Diminishing Lax-Fiedrichs 
(TVDLF) scheme [H] for relativistic MHD applications. Temporal second 
order accuracy is achieved by the Hancock predictor step 



..n+l/2 ^ jrn 1 ^t 

' ' 2Aa; 



Fiur + -Af/n - F(f/r - -mi: 



(29) 



Here, AU denote the cell-center to cell-face limited slope used in the TVDLF 
scheme. In our work we will mostly use the rather diffusive, but stable 'min- 
mod' limiter 



AUi = sgn{Ui - Ui_i) max [0, min {| Ui - Ui^i |, 

(f/,+i-f/i)sgn([/,-[/,_i)}]. (30) 

In the full correction step, the numerical fluxes are 



— \ r I = — 



^."h-'^'H 



(31) 



where the left and right states are 



u:i.=ur:,'^'-Au::,'^'/2, (32) 

respectively. This TVDLF scheme does not use a Riemann solver. The only in 



formation needed is the fastest characteristic wave speed Cmax = niax( 



AF 



A^ 



In this second order scheme, some improvement is obtained if the limited 
slopes are calculated via the primitive variables. We have best experience with 
(p, rv,p, B). Note the inclusion of the Lorentz factor F in the fluid velocity. 
We also experimented with HLL and HLLC solvers (see e.g. [32]), which use 
more information of the Riemann fan at cell interfaces, but leave their ap- 
plication outside the scope of this paper. It is our impression that the use of 
grid adaptivity makes the difference between these base solvers of secondary 
importance. 



3.4 Parabolic magnetic monopole treatment 



In our multidimensional RMHD applications, we handle the V ■ B = con- 
straint by adding a diffusive source term proportional to V V • B to the induc- 
tion equation. The diffusion coefficient is determined by setting the maximum 
allowed diffusion time step equal to the CFL time step. On a cartesian mesh, 
we obtain the source term update 




(33) 

where < Cj < 2 for stability reasons. For most applications, a value Cd = 1 
is advised. This is similar to the strategy used for non-relativistic MHD on 
curvilinear grids in [^ . 



This type of treatment for the V ■ B = can be regarded as the parabolic vari- 
ant of the hyperbolic/parabolic treatment discussed for ideal, non-relativistic 
MHD in Dedner et al. [10] . Our source treatment redistributes potential local 
monopole errors over a wider area than where they would normally concen- 
trate. It should be noted that this is more meant as a means to stabilize the 
overall numerical scheme and to avoid potential numerical instabilities, than 
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as a manner to annihilate discrete monopole contributions alltogether. The 
latter is not needed in any numerical integration, where we will always have 
truncation errors in the magnetic field components, even if we numerically 
ensure a kind of minimal V ■ B. The hyperbolic cleaning from [10] in addition 
advects these locally occuring numerical monopole errors, while damping them 
in a similar fashion. In case of shock-dominated problems (like those presented 
further on), discrete monopole errors are continuously arising at shock loca- 
tions, so the damping by itself is in our opinion the most important part of the 
error control strategy. The same idea was first suggested in electromagnetic 
PIC codes, by Marder [2H], and is now also routinely used [39] in the only 3D 
global tokamak plasma simulations feasible to date, performed by the NIM- 
ROD consortium. One scheme for resistive relativistic MHD has recently been 
presented in [2S], where the fiuxes are handled using an HLL scheme instead 
of our TVDLF method, while the divergence treatment there uses the Gener- 
alized Lagrange Multiplier approach from Munz et al. [M]. This in essence is 
the hyperbolic variant from [10], and our parabolic treatment can be seen as 
belonging to the same family of monopole treatments. In [22] , a comparison of 
the parabolic treatment against other popular source term strategies (Powell 
source terms as e.g. described in [36j, or only modifying the induction equation 
as suggested by Janhunen [21]) was performed for various multidimensional, 
non-relativistic MHD problems using AMR. The most extensive comparison of 
magnetic divergence control in mult i- dimensional, non-relativistic MHD sce- 
narios for fixed grids can be found in Toth j40j . In particular, the Powell source 
term treatment was then tested against a variety of constrained transport im- 
plementations, which insist on ensuring V ■ B = to machine precision in one 
particular pre-chosen discretization only, as well as against an elliptic clean- 
ing scheme. The latter strategy projects the obtained B* to a divergence free 
magnetic field B = B* — V0, and involves the solution of a Poisson equation 
for (p. Noteworthy is that it was also found in [ID] that one does not need 
to enforce its solution to machine precision either. All of these schemes were 
found to yield acceptable simulation results on the nine tests verified there. 
Some of these tests were revisited in the grid adaptive computations presented 
in [22j , where only different source term treatments were compared. Those that 
specifically only modify the induction equation are readily incorporated in the 
relativistic MHD regime, and they do not violate the conservation of the other 
than magnetic variables. In our multidimensional tests below, we quantify the 
remaining errors in V ■ B for illustration purposes only: the fact thay they 
remain bounded at all times is the most important observation to be checked 
there. 
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4 Adaptive Mesh Refinement 



There are two AMR versions implemented in our AMRVAC code, namely a 
modified version of the original patch approach of Berger [5] and a hybrid 
block-based j32| strategy. We will give a brief outline of these methods and 
refer to the literature for details. 

Once a procedure is given for detecting cells that are needed for resolving 
fiow features, the AMR must arrange these cells in a hierarchy of properly 
nested grids. In the original patch-based scheme of Berger, each cell fiagged 
for refinement is surrounded by a buffer- zone of a user given size to ensure that 
discontinuities and other regions that need high resolution do not propagate 
to coarser cells. This collection of flagged cells are then changed so that it 
satisfles the proper nesting: each level / cell is in between level I + 1 and 
I — 1. These cells are subsequently stored into a hierarchy of properly nested 
grids (patches). The patches are then bisected till a given efficiency is reached. 
This efficiency is expressed as the ratio of the fiagged to the total amount of 
cells within a patch. Finally, a patch merging process is called to reduce the 
computational cost of too many small grids. In the AMRVAC code, see |22] , 
the overlap of patches on the same AMR level is avoided, while a minimal 
efficiency is enforced. 

In the same AMRVAC code, yet another AMR scheme is implemented. This 
so-called hybrid block-AMR method |12] uses, like block-based AMR tech- 
niques, an equal number of grid cells for each grid in the entire grid hierarchy. 
The basic structure is an octree (in three dimensions). However, if a block is 
fiagged for refinement, this scheme relaxes on the standard approach where 
a block triggers in a D-dimensional calculation 2^ new sub-blocks (children). 
This consequently introduces incomplete block families in the grid hierarchy. 
Therefore, the hybrid method approaches the optimal fit of the grid structure 
in the patch scheme. However, due to fixing the number of cells per grid, the 
good cache performance of the common tree block-based approach is fully 
exploited. 

The procedure to identify which cells are to be triggered for refinement relies 
on a type of Richardson extrapolation. The error estimator implemented in 
AMRVAC is a variant of the procedure given in [6]. Given the solution vector 
f/"~^ and up on level / and with a time difference At"~^, the error estimator 
will first determine ?7;"lV i^ ^^^ following two ways: 

• Coarsen [/""-^ and then advance to time t = t"^-^ + 2At""\ 

• Advance Up to time t = t"""*^ + 2At"~^ and then coarsen. 

The resulting solution vectors are compared and, based on a certain selection of 
the conservative variables with mutual weighting factors, cells will be fiagged 
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in an automated fashion for refinement if a given tolerance is exceeded. In 
addition, user-enforced refinement is possible and the auxiliary variables, like 
the Lorentz factor, can be exploited in the error estimator. 

Essential in the AMR strategy is the restoring of the global conservation across 
the entire grid hierarchy. This amounts to fixing the fluxes of coarse cells with 
the fluxes of the neighboring fine cells. Moreover, the prolongation, restric- 
tion, and/or temporal interpolation between different grid levels need to be 
performed on the conservative variables. To avoid the possible introduction of 
negative pressure, especially in the vicinity of small ratio of kinetic to mag- 
netic pressure, we have implemented in AMRVAC the option to use primitive 
variables during the regridding process. Another option, as pointed out earlier, 
would be to switch to the entropy DS instead of energy during the regridding 
process. Then the AMR scheme is still conservative, but avoids the introduc- 
tion of negative pressure. 



5 Numerical results 



5.1 Riemann problems 



A direct validation of the code is provided by performing a series of Riemann 
problem tests, where an initial discontinuity is left to evolve dynamically. For 
relativistic MHD, a recent contribution by [17] documented how one can ob- 
tain the exact 'analytic' solution when accounting for the up to 7 wave signals 
that may emerge out of the t = problem. The method was demonstrated 
for 10 specific initial conditions, and we reproduce all these 10 cases numeri- 
cally here. They collect various tests reported in recent code developments for 
relativistic MHD p^lH] . augmented with some new Riemann problems. Our 
results are shown in Figs. [l]J2| where we overplot in all cases the exact solution 
generated by [TTj. We briefly comment on each case in what follows. We typ- 
ically compute with up to 8 grid levels, and use a base resolution of 60 unless 
stated otherwise. The exact initial conditions are fully specified in [17], and 
we refer to that paper for details. 

The first two tests have 3^ = throughout, so that only two fast signals and 
a tangential discontinuity emerges. The first two rows of Fig. [T] correspond 
to the Shock- Tube 2, and Generic shock tube test, respectively. The first has 
a left-going fast rarefaction, a tangential discontinuity, and a fast shock, and 
presents no major difficulty. The generic shock tube test has a leftgoing fast 
shock and rightgoing rarefaction, with a tangential discontinuity in between. 
We find overshoot errors at both the shock and the tangential discontinuity, 
which diminish only by raising the overall resolution significantly: in Fig.jTjthis 
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test here exploited 240 base level grid points, and we still have an erroneous 
variation in between shock and tangential signal in Vz- 

The next set of 8 tests considers cases where Bx ^ Q. The third to fifth row 
in Fig. [l] shows the outcome for the 3 coplanar problems, where at most 5 
wave signals can emerge. The third row has no y— or z— velocity nor magnetic 
field components, and represents a typical challenge when a contact discon- 
tinuity is in close vicinity to a fast shock. Our AMR computation with the 
TVDLF scheme adequately resolves the narrow structure, with a rather un- 
avoidable large number of grid points to represent the contact. The fourth 
test is a collision leading to two pairs of left and rightgoing fast and slow 
shocks. Similar to all documented numerical solutions, we can hardly avoid 
to generate a central error in density p, which should remain constant in be- 
tween the two slow shocks. The final case shown in Fig. [T] is the relativistic 
analogue of the Brio- Wu test [8] . A left going rarefaction, a slow compound, 
a contact discontinuity, and a slow and fast shock are encountered from left 
to right. The correspondence with the analytic solution is satisfactory, with 
many grid points representing the contact jump. Note that the method to 
generate the 'exact solution' excludes the possibihty of compound waves, in 
part invalidating the comparison there. 



The tests collected in Fig. p^ reproduce tests from [1] and the generic Alfven 
test introduced by [T7] . The top row has left going fast and slow rarefactions, 
a contact and rightgoing slow and fast shocks. Only the many grid points 
in the contact is arguable, but at the same time the variation is captured 
at the highest grid level activated. The next test has a similarly structured 
outcome, with an even more extreme length contraction effect at play between 
right-going contact, slow and fast shock. The latter are accurately captured 
at the correct amplitude, as best seen in the By plot. The third row of Fig. [2] 
yields another stringent collision test analogous to the one shown in Fig. [ij 
With no major differences in solution strategy, the numerical result here (with 
base resolution 120) is far less polluted by a central density error. The last 
two tests trigger all 7 wave signals from the initial discontinuity, with Alfven 
discontinuities in between the slow and fast signal pairs. At the times shown, 
the spacing between Alfvenic and slow signals can be very close still. For the 
fourth row in Fig. |2j the Balsara test 5 case, it is expected to be down to 
order 0.001 for the leftgoing Alfven discontinuity and slow rarefaction. It is 
seen in the Bz plot how even higher effective resolution would need to be used 
to capture this transition exactly. Still, we correctly find all wave signals. The 
last generic Alfven test has a similar challenge for the Alfven signals adjacent 
to slow shocks. The outmost left-going fast rarefaction and right-going fast 
shock are captured on the lowest grid resolution only for this test, which 
can be changed by fine-tuning the employed refining criterion. In summary, 
our grid-adaptive numerical solutions show a favorable agreement with the 
analytic results in all cases. 



14 



5.2 Multidimensional tests 



In a first 2D test, we recompute the relativistic rotor problem, a test first 
performed in non-relativistic MHD settings, and subsequently modeled in a 
relativistic variant by [12] . We simulate the very same problem on a larger do- 
main [0,2]^, in order to follow the rotor evolution at higher effective resolution 
to a longer time than presented originally. The polytropic index is constant at 
7 = 5/3, while initial pressure is unity throughout. A high density circular disk 
with p = 10 rotates anti-clockwise at uniform angular velocity u = 9.95 within 
radius 0.1 from the center of the domain. The disk discontinuously connects 
to a tenfold lower density, static medium. The entire domain is pervaded by a 
homogeneous horizontal field B = e^;. Using 7 refinement levels we acquire an 
effective resolution of 6400 x 6400. This represents locally a four times higher 
resolution than used in [12j, which compensates for the difference in order of 
accuracy employed (a third order method versus our second order scheme, on 
smooth solutions). 

We intentionally followed the rotor evolution to twice the time reported earlier, 
till t = 0.8. In [12], slight corrugations in density could be detected cospatial 
with shear flow regions at time t = 0.4, and we intended to investigate their 
potential role in any further nonlinear evolution. Snapshots at t = 0.4 and 
t = 0.8 are shown in Fig. [3j Fast shock fronts can be detected to travel out- 
wards into the static surroundings, and inwards towards the disk center. Slow 
rarefactions are found in between, and the fleld deflections in effect brake the 
initial fast rotation (F ^ 10). The overall evolution remains pointsymmetric. 
At t = 0.4 the fleld in the disk has rotated over about 90°, and the in- and 
outwards traveling shocks can be clearly detected. We found somewhat higher 
instantaneous Lorentz factors than those reported before, and no evidence of 
density corrugations. Also, as shown in the second snapshot where the in- 
ward shocks have already collided, no clear indication of a shear-induced flne 
structure was found. The automated grid reflnement does follow the density 
variations at the highest grid resolution throughout the computation. An im- 
pression of the location of the intermediate level 5 grids (two more grid levels 
exist on top of this level) is shown in Fig. p] 

To illustrate the magnetic monopole control used in this particular simulation, 
we provide various quantitative measures for it in Fig. |4j As stated earlier, the 
parabolic treatment is meant to stabilize the computation and uses a discrete 
evaluation of the magnetic fleld divergence in a diffusive type source term. The 
errors themselves are unavoidably created continuously at the location of the 
strongest discontinuities. In this simulation, even the flrst timestep introduces 
flnite monopole errors at the border of the 'rotating disk', since we effectively 
break the fleld lines there (e.g. the central fleld line gets disconnected in two lo- 
cations, where from one grid point to the next we jump from static to Lorentz 
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factor 10 flow). These inonopoles can locally and temporarily be of order unity, 
while the diffusive treatment then spreads and diffuses these errors during the 
computation. In Fig. [4], the left panel shows two domain averaged error mon- 
itors as a function of time: it can be seen that these mean values remain at 
(9(0.01) in this computation. The fact that they do not grow without bound 
is the most important observation, conflrming their stabilizng role. The right 
panel from Fig. |4]is at the flnal time t ~ 0.8, and shows the instantaneous dis- 
tribution of the monopoles, as evaluated in centered difference approximation. 
Black values are locally order unity, and the largest errors necessarily coincide 
with the various shock wave fronts. Note that our restriction and prolongation 
strategies do not particularly enforce solenoidal fields in any discrete sense, so 
grid level boundaries can temporarily be detected in such error maps. Again, 
the role of the parabolic term then acts to diffuse such errors at their maximal 
rate and it is important to note that this does not affect any of the employed 
conservative variables adversely. 

In a second 2D application, we generalize a shock-dominated time-dependent 
problem frequently used in benchmarking classical MHD codes to the rela- 
tivistic regime. The non-relativistic test [35] considers a Mach 1 vortex super- 
posed on a multiple magnetic island configuration, on a doubly periodic [0, 2tt]'^ 
Cartesian domain. It starts from uniform density and pressure throughout, and 
the supersonic rotation concentrates magnetic field gradients in thin, localized 
current sheets from which shock fronts originate, which subsequently interact. 
Our relativistic analogue considers a relativistically hot gas, where the internal 
energy dominates over the rest mass contribution, such that the relativistic 
sound speed approaches its maximal value Cg ~ y/j — 1. This is consistent 
with a polytropic index value 7 = 4/3, and we set the initial pressure p = 10 
while proper density p = 1. The vortex imposes a velocity field 

\ = —Asm{y)ex + Asm{x)ey, (34) 

where A = 0.99 /a/2, ensuring subluminal and supersonic velocities. The maxi- 
mal initial Lorentz factor is then about 7. The magnetic field is then initialized 

at 

B = — sin {y)ex + sin {2x)ey, (35) 



which makes the ratio of magnetic to thermal pressure attain a maximal value 
of (3i = 0.098. We simulated this problem for times beyond t = 12.5, at which 
time the maximal Lorentz factor encountered has dropped to about F = 1.526. 
The simulation used a 40 x 40 base resolution, with a maximal 7 refinement 
levels, effectively mimicking a 2560 x 2560 resolution. 

The initial magnetic topology is characterized by alternating X and O nul- 
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points (where the field vanishes), with 4 different X and 4 O type nuls in 
the doubly periodic domain. Along the y = vr horizontal, we thus encounter 
two islands of closed field (to the left and right of the central X point), and 
this pattern repeats on y = 0, 27r with distinct 7r/2 phase difference. The su- 
perposed vortical velocity will immediately displace the left-central island in a 
diagonally upwards fashion, while distorting the right-central island diagonally 
downwards. The double periodicity implies that the island situated midway 
the horizontal boundaries gets squashed in the process, and oriented in a diag- 
onally downwards current sheet. This violent compression will drive two shock 
waves, one from either side of this sheet, traveling against the original flow 
direction. In a similar fashion, the flow induces a strong shearing at the X 
point midway the vertical sides, with similar accumulation of matter along 
a diagonally upward pointing sheet. There too, two shock fronts separate off 
the compression zone. The flrst panel of Fig. [5] at t = 2.82 superposes the 
field structure on the (logarithmic) proper density, clearly showing the four 
diagonal shock fronts and the island deformations just described. These four 
interacting shock fronts meet up in the centre of the domain, while the sheet 
formed in the compressed O configuration eventually demonstrates a spon- 
taneous break-up forming a series of islands. This tearing type reconnection 
happens at about t = 4.6. In this sudden topological magnetic reconfiguration, 
some of the smaller islands get accelerated towards the shock fronts traveling 
away from the sheared and compressed X point. In the second panel of Fig. [5] 
at t = 6.85, some of these islands and the shock front deformations they cause 
can be seen. These and similar interactions occuring with the diagonal shock 
fronts converging to the sheared X sheet, drive a second sequence of strong 
density variations colliding upon both the compressed O point (now with its 
island structure) and the sheared X sheet. A second quadruplet of shocks 
then propagates away from these locations. This causes intricate shock-shock 
interference patterns with the original 4 shock fronts. The sheared X sheet 
in fact also gets torn up into an island chain structure, as seen in the third 
panel from Fig. |5| Eventually, also this second shock sequence meets up at 
the centre, while only some of the islands from the original O and X sheet 
tearing events survive as localized density enhancements or depletions. The 
final panel of Fig. |5] shows the by now rather complicated density distribu- 
tion near time t = 12.5. The point-symmetry about the centre of the domain 
is preserved perfectly throughout the entire grid-adaptive computation. Note 
that this long-term computation gives strong evidence for intricate reconnec- 
tion events, which will be infiuenced by the numerical scheme employed: it 
could be of interest to benchmark several higher order schemes, in combina- 
tion with varying strategies for magnetic monopole control, on this problem 
in particular. Ultimately, deviations from perfect conductivity would need to 
be explicitly accounted for. 

As a final illustration of the parabolic monopole treatment used in this test. 
Fig. |6] collects the temporal evolution of various domain averaged error moni- 
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tors over most of the simulated period. Once more, the error diffusion approach 
works as intended, despite the presence of very strong interacting shock fronts, 
the sudden appearance of magnetic island chains, and the continuously adjust- 
ing AMR grid hierarchy. 



6 Conclusion 



We provided details on our relativistic grid-adaptive MHD code, and tested it 
using recently available Riemann problem solutions, as well as multi-dimensional 
setups. The latter include a new long-term simulation demonstrating shock- 
shock interactions and reconnection events, and this could be of interest to 
benchmark existing shock-capturing algorithms on long-term shock-dominated 
relativistic magnetized flow problems, such as those recently presented in |25] . 

In future work, the RMHD code AMRVAC will be used to model AGN jet 
propagation [23] and GRB dynamics. The magnetic field is suggested as possi- 
ble ingredient to achieve the high Lorentz factor flows, such as reached by the 
GRB fireball, and plays a crucial role in both AGN and GRB flow collimation. 
Other future projects concentrate on implementing a more realistic equation 
of state to investigate the launching and the propagation of relativistic jet. 
We will discuss these astrophysical applications with more detail in future 
publications. 
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Fig. 1. RMHD Riemann problems. Crosses indicate the AMR solutions, while solid 
lines are the exact solutions. 
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Fig. 2. RMHD Riemann problems. 




Fig. 3. The relativistic rotor at times t = 0.4 and t = 0.8, showing at left the Lorentz 
factor and the location of the intermediate level 5 grid blocks, and the magnetic 
field lines on top of proper density distribution at right. 
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